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Numerical  simulation  of  protoplanet  ary  vortices 

By  H.  Lin,  J.A.  Barranco  f and  P.S.  Marcus  $ 


1.  Introduction 

The  fluid  dynamics  within  a protoplanetary  disk  has  been  attracting  the  attention  of 
many  researchers  for  a few  decades.  Previous  works  include,  to  list  only  a few  among 
many  others,  the  well-known  a-prescription  of  Shakura  & Sunyaev  (1973),  the  convective 
and  instability  study  of  Stone  &;  Balbus  (1996)  and  Hawley  et  al.  (1999),  the  Rossby  wave 
approach  of  Lovelace  et  al.  (1999),  as  well  as  a recent  work  by  Klahr  & Bodenheimer 
(2003),  which  attempted  to  identify  turbulent  flow  within  the  disk.  The  disk  is  commonly 
understood  to  be  a thin  gas  disk  rotating  around  a central  star  with  differential  rota- 
tion (the  Keplerian  velocity),  and  the  central  quest  remains  as  how  the  flow  behavior 
deviates  (albeit  by  a small  amount)  from  a strong  balance  established  between  gravi- 
tational and  centrifugal  forces,  transfers  mass  and  momentum  inward,  and  eventually 
forms  planetesimals  and  planets. 

In  earlier  works  (Barranco  & Marcus  2000;  Barranco  et  al.  2000;  Lin  et  al.  2000)  we 
have  briefly  described  the  possible  physical  processes  involved  in  the  disk;  we  have  pro- 
posed the  existence  of  long-lasting,  coherent  vortices  as  an  efficient  agent  for  mass  and 
momentum  transport.  In  particular,  Barranco  et  al.  (2000)  provided  a general  mathe- 
matical framework  that  is  suitable  for  the  asymptotic  regime  of  the  disk;  Barranco  & 
Marcus  (2000)  addressed  a proposed  vortex-dust  interaction  mechanism  which  might 
lead  to  planetesimal  formation;  and  Lin  et  al.  (2002),  as  inspired  by  general  geophysi- 
cal vortex  dynamics,  proposed  basic  mechanisms  by  which  vortices  can  transport  mass 
and  angular  momentum.  The  current  work  follows  up  on  our  previous  effort.  We  shall 
focus  on  the  detailed  numerical  implementation  of  our  problem.  We  have  developed  a 
parallel,  pseudo-spectral  code  to  simulate  the  full  three-dimensional  vortex  dynamics  in 
a stably-stratified,  differentially  rotating  frame,  which  represents  the  environment  of  the 
disk.  Our  simulation  is  validated  with  full  diagnostics  and  comparisons,  and  we  present 
our  results  on  a family  of  three-dimensional,  coherent  equilibrium  vortices. 


2.  Governing  equations 

Scaling  issues  have  been  discussed  in  both  Barranco  et  al.  (2000)  and  Lin  et  al.  (2002), 
and  shall  not  be  repeated.  Here  we  simply  list  the  governing  equations  that  we  propose 
to  solve  numerically,  namely 

+ p(v  • Vv)  = -2pQ.0z  x (v  - v)  — Vp  - pft20zz,  (2.1) 


dT 

— + v • VT  = 


V • (pv)  = 0, 


Tradfy) 
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(2.2) 

(2.3) 
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Here  the  equations  are  written  in  a rotating  frame  with  an  angular  velocity  fi0,  x is  the 
streamwise  ( <f> ) direction,  y is  the  radial  (r)  direction,  z is  the  direction  about  which  the 
frame  is  rotating,  and  z is  its  corresponding  unit  vector.  In  our  notation  system  p , p,  T0 
and  v denote  the  base  state  density,  pressure,  temperature  and  velocity,  respectively,  and 
their  forms  will  be  given  shortly.  The  tilde  quantities  p and  p denote  the  deviations  of 
density  and  pressure  from  the  base  states,  respectively,  and  7 is  the  ratio  of  specific  heats. 
Finally  H0  is  the  scale  hight  of  the  disk  and  Traj  is  a radiation  time  scale.  Compared 
with  the  governing  equations  (4. 1-4.4)  of  Lin  et  al.  (2002),  we  have  further  expanded 
the  energy  equation  making  use  of  the  anelastic  continuity  equation  (2.2);  we  have  also 
linearized  the  equation  of  state  to  obtain  equation  (2.4)  which  is  consistent  with  our 
current  low-Mach  number  asymptotic  regime. 

We  nondimensionalize  equations  (2. 1-2.4)  with  the  following  scales: 

1 £]  = »«  M = ‘"  ['!  = $ = £■ 


[p]  = Po,  [p]  = PoC2s  = Po , [T}  = To, 

where  cs  is  the  speed  of  sound,  p0  is  a characteristic  density,  and  T0  is  the  characteristic 
(base)  temperature.  The  resulted  dimensionless  equations  are 


= -p(v  • V)v  - 2pz  x (v  - v)  - Vp  - pzz, 


dT 

dt 


V • (pv)  = 0, 

= —v  • VT  — (7  — 1)  (l  + ^) 


VzZ  - 


tVad(^) 


p = p + pf. 


(2.5) 

(2.6) 

(2.7) 

(2.8) 


The  base  states,  as  discussed  in  Lin  et  al.  (2002),  are  given  (in  dimensionless  form)  as: 


V = P = Po(y)exp 


(4). 


(2.9) 


3 

v = -yx,  (2.10) 

where  pc(y ) represents  the  radial  variation  of  the  base  density,  and  x is  the  unit  vector  in 
the  x direction.  The  dimensionless  base  temperature,  which  is  assumed  to  be  constant, 
simply  becomes  Ta  = 1.  For  the  radiation  time  rrad , we  currently  assume  the  following 
dimensionless  form, 

rTad{z)  = r0exp  (-z2)  . 

We  set  t0  to  be  approximately  106  times  that  of  the  orbital  period  (denoted  rorb).  When 
such  a value  is  assumed  there  is  almost  no  radiative  dissipation  at  the  midplane,  and 
the  radiation  time  scale  is  comparable  to  the  orbital  period  close  to  the  upper  and  lower 
boundaries  of  the  disk.  The  radiation  scheme  captures  only  the  long  term  dissipative 
behavior  of  the  disk,  and  has  a negligible  effect  on  dynamics  that  have  time  scales  com- 
parable with  orbital  periods. 

Equations  (2. 5-2.8),  together  with  appropriate  boundary  conditions,  are  the  set  of 
equations  we  solve  numerically. 
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3.  Numerical  scheme 

Before  presenting  in  detail  the  numerical  algorithm,  it  is  appropriate  to  address  a few 
of  the  special  challenges  involved  in  the  system  we  are  trying  to  solve.  First,  we  have 
a base  shear  flow  in  a rotating  frame  as  indicated  by  equation  (2.10).  For  this  kind  of 
problem  it  is  natural  to  think  about  a sliding-box  formulation,  as  it  has  been  implemented 
by  Rogallo  (1981)  and  Hawley  et  al.  (1999).  Two  of  the  current  authors  (J.A.B.  and 
P.S.M.)  have  also  developed  a similar  algorithm  to  solve  the  system  of  equations  (2.5- 
2.8)  using  the  sliding-box  methodology.  In  this  sliding-box  formulation,  flow  in  the  radial 
direction  is  artificially  periodic  after  subtraction  of  the  convection  ensued  from  the  base 
shear  flow,  and  the  vertical  direction  is  represented  in  Chebyshev  series  to  capture  the 
non-periodic  buoyancy  (gravity)  force.  Nonetheless,  because  in  the  current  code  we  aim 
to  represent  a non-periodic  radial  thermodynamic  (potential-vorticity)  background  (see 
discussions  in  Lin  et  al.  (2002)),  we  can  no  longer  assume  the  artificial  periodicity  in  y, 
and  functions  in  this  direction  should  be  instead  represented  in  Chebyshev  series.  This 
complication  forces  us  to  assume  an  artificial  periodicity  in  the  vertical  direction  under 
the  current  framework  of  a Fourier-Fourier-Chebyshev  representation,  and  we  shall  have 
more  discussions  regarding  this  point  later.  Secondly,  we  account  for  density  variation 
(albeit  with  a simplified  form)  via  the  anelastic  continuity  equation  (2.6),  and  capture  the 
stably-stratified  physics  through  a vertically  decaying  density  p,  a temperature  (energy) 
evolution  equation  (2.7),  and  a buoyancy  (gravity)  force  term  in  the  momentum  equation 
(2.5).  Thirdly,  because  the  base  density  has  the  fast-decaying  form  of  equation  (2.9)  ( e.g . 
at  z = ±3 H0,  the  base  density  decays  to  only  1%  of  its  value  at  the  midplane),  it 
should  be  carefully  treated  to  avoid  numerical  instability.  Finally,  the  coherent  vortices 
we  are  trying  to  identify  assumes  full  three-dimensional  structures  for  reasons  discussed 
in  Barranco  et  al.  (2000). 

We  solve  equations  (2.5-2.8)  with  a pseudo-spectral  scheme.  We  represent  functions  in 
the  x and  z direction  with  Fourier  series,  and  functions  in  the  y direction  by  Chebyshev 
series.  Because  of  the  artificial  periodicity  in  the  z direction,  we  have  to  replace  the 
variable  z in  the  buoyancy  force  term  pz  with  a new  periodic  function  z*  that  equals  z 
for  most  of  the  domain,  but  comes  back  to  zero  to  enforce  periodicity  toward  the  upper 
and  lower  boundaries  (see  Figure  1).  Because  the  base  density  p (and  hence  momentum, 
energy,  etc.)  becomes  very  small  where  z*  has  the  largest  deviation  from  z,  and  because 
our  vortices  have  a compact  structure  around  z = 0,  the  influence  of  this  replacement  of 
z with  z*  is  assumed  to  be  minimal. 

For  time  integration  we  adopt  a Leap-Frog  scheme  for  calculation  of  the  nonlinear 
terms.  Our  algorithm  can  be  compactly  written  as 

uM+ 1 = uM-1  + 2At(nM  - pMz* z),  (3.1) 

vv^v.  *M+1.  <3-2> 

UM+1  _ £M+ i _ 2A tVpM,  (3.3) 

tm+i  = (Tm- i + 2A  tnT)e-2At/T™*,  (3.4) 

where 

u = pv  (3.5) 

is  momentum  and  our  working  variable,  the  superscript  M denotes  quantities  at  the  Mth 
time  step,  n is  the  sum  of  the  nonlinear  terms  on  the  RHS  of  equation  (2.5)  excluding  the 
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Figure  1.  The  artificially  enforced  periodic  function  z * (solid),  compared  with  the  original 
linear  function  z (dash-dotted).  In  this  plot  the  vertical  domain  extends  from  —3 H0  to  3 H0\  the 
value  of  the  base  density  p at  the  upper  and  lower  boundaries  (z  « ±3 Ha)  is  about  1%  of  its 
value  at  the  midplane  [z  = 0).  Our  vortex  structure,  as  we  shall  present  later,  is  compact  in  z 
and  primarily  assumes  the  domain  from  —H0  to  Ha. 


buoyancy  force,  and  nT  is  the  sum  of  the  nonlinear  terms  on  the  RHS  of  equation  (2.7) 
excluding  the  radiation  term.  The  Helmholtz  equation  (3.2)  is  inverted  to  solve  for  pres- 
sure which  consequently  enforces  the  divergence  free  condition  (2.6)  at  the  ( M+l)th  step. 
In  our  code  we  have  also  implemented  a r-method  that  ensures  equation  (2.6)  is  satisfied 
for  all  Chebyshev  modes.  We  found  that  the  latter  method  improves  numerical  stability. 
For  boundary  conditions,  we  are  currently  implementing  no-flow  condition  at  y = ±1 
(which  is  enforced  by  the  pressure  boundary  condition  when  solving  for  pM  in  (3.2)). 
The  advantage  of  doing  this  is  to  contain  our  flow  and  enforce  conservation  laws  within 
the  computational  domain.  Finally  we  apply  a 6th  order  hyperviscosity/hyperdiffusion 
on  both  the  momentum  and  energy  equations  to  ensure  numerical  stability. 

Note  that  we  have  not  adopted  the  more  commonly  practiced  Adams-Bashforth  (for 
nonlinear  terms)  and  Crank-Nicolson  (for  pressure)  time-integration  scheme.  This  is  be- 
cause the  current  Leap-Frog  scheme  has  better  performance  in  terms  of  avoiding  the 
instability  caused  by  the  large  density  variation  of  p in  the  vertical  direction.  Note  also 
that  when  inverting  the  Helmholtz  problem  (3.2)  for  pM , the  uM+1  term  on  the  RHS 
contains  a buoyancy  term  pM  z*,  in  which  the  density  pM  is  related  to  the  pressure  pM 
via  the  equation  of  state  (2.8).  This  coupling  problem  is  resolved  by  an  iterative  method. 
We  use  pM~l  as  the  initial  value  pft1  to  enter  the  iteration  on  the  RHS  of  equation  (3.2). 
Once  the  pressure  p is  obtained  by  inverting  the  Helmholtz  problem,  it  is  used  to  obtain 
a new  density  function  using 

p™x=p¥-pTM,  (3.6) 

which  subsequently  enters  the  Helmholtz  solver  as  the  new  density,  we  found  that  ~ 5 
times  of  this  procedure  give  excellent  convergence  at  an  acceptable  computational  cost. 

To  improve  the  speed  of  the  code  we  have  also  implemented  a “faster”  variation  of  our 
algorithm.  This  implementation  is  based  on  the  following  consideration.  We  can  rewrite 
the  momentum  and  energy  equations  as 

du  du  _ „ 


(3.7) 
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dT  dT  T 

- = -^+„r-  — 

where  a = | is  the  shear  rate  of  the  base  flow,  and  n and  are  the  same  as  previ- 
ously defined.  In  this  case,  however,  we  further  exclude  the  base-shear  convection  terms. 
The  base  shear  flow  poses  a major  constraint  on  the  size  of  the  time-step  if  calculated 
explicitly,  whereas  the  flow  structure  we  are  interested  in  is  only  a small  deviation  from 
this  base  flow  (Barranco  et  al.  2000).  We  therefore  separate  the  base  convective  terms 
— o-ydn  anci  —uyQI-  and  treat  them  analytically.  Our  new  algorithm,  in  which  the  base 
convection  is  represented  in  terms  of  waves,  is  written 


uM+1  = e~2icAt UM-'  + e~icAt2At(nM  - VpM_1  - pM~1z* z), 


(3.9) 


UM+1  _ flM+i  _ 2A t[V(pM  - PM_1)  4-  ( PM  - pM~1)z*z),  (3.10) 

TM+1  = (e-2ic*tTM-l  + e-icAt2£tnM}e-2At/Trad ' (3.H) 

Here  c = kxay,  that  is,  these  exponents  are  applied  for  each  wave  number  kx  in  the 
x direction  and  at  each  point  y in  the  transverse  domain.  This  formula  assumes  the 
form  of  a semi-implicit,  predictor-corrector  algorithm  of  second-order  accuracy.  Again  a 
Helmholtz  problem  is  solved  to  obtain  pressure  and  enforce  continuity  equation,  and  an 
iterative  method  is  used  to  simultaneous  obtain  pM  and  pM . In  this  new  algorithm  our 
effective  flow  number  is  calculated  from 


|v|  = |v- v|  < |v|. 

Because  the  maximum  value  of  v is  typically  only  1/5  that  of  v,  we  can  use  significantly 
larger  time  steps. 

Our  algorithms  are  implemented  for  parallel  computation  using  MPI,  with  a com- 
munication algorithm  for  FFT  provided  by  Dr.  Alan  Wray.  We  test  and  run  our  code 
on  the  SGI  Origin  2000  clusters  at  NASA  Ames.  We  have  performed  diagnostics  such 
as  conservation  of  mass,  momentum  and  energy.  Aside  from  the  results  that  we  shall 
present  below,  we  have  also  tested  finer  resolutions  for  convergence.  We  have  followed 
the  standard  diagnostic  procedures  and  the  details  shall  not  be  presented  here. 


4.  Results 


We  now  present  our  computational  results  for  three  dimensional  vortices.  We  shall 
be  seeking  the  existence  of  three-dimensional,  coherent  vortices  which  can  survive  in  an 
environment  such  as  that  of  the  protoplanetary  disk.  We  will  call  these  vortices  henceforth 
“equilibrium”  vortices. 

Because  we  seek  vortices  that  are  in  dynamic  equilibrium  with  their  environment, 
a natural  approach  would  be  to  generalize  the  idea  of  the  well-known  Moore-Saffman 
vortices  (see,  e.g.  Saffman  (1992))  to  three  dimensions.  In  Lin  et  al.  (2002)  we  have 
discussed  columnar  Moore-Saffman  vortices,  which  are  quasi-two-dimensional  objects.  In 
the  following  we  present  results  for  three-dimensional  vortices,  which  resemble  Moore- 
Saffman  vortices  in  the  horizontal  planes,  but  have  compact  support  in  the  vertical 
direction.  Before  we  start  it  is  useful  to  list  the  vorticity-aspect-ratio  relation  given  by 
Moore-Saffman  theory: 


= A 


A + 1 


(4.1) 


(7  A — 1 ’ 

where  u>z  is  the  strength  of  the  constant  vorticity  patch  superimposed  on  the  shear,  and 
A < 1 is  the  ratio  of  the  semi-axes  of  the  elliptic  patch. 
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Figure  2.  Three-dimensional  vortex  structure  in  the  absence  of  shear.  Here  we  show  the  iso- 
surface at  the  value  u>z  = —1,  and  at  the  time  t = 20  rort,.  Our  computational  resolution  is 
128  x 64  x 96  on  a computational  domain  of  8 x 2 x 6 (not  entirely  shown  in  this  graph). 


4.1.  Vortices  without  shear 


We  first  present  a relatively  simple  case  in  which  we  set  the  shear  to  be  zero,  that  is, 
a = 0.  This  case  will  facilitate  a comparison  with  our  later  results  where  we  see  the 
effects  of  a strong  shear  on  the  vortex.  As  the  initial  condition,  we  specify  a vorticity 
function  as 


& z 


w0exp(--y)’  yj x2  + y2  < r0 
0,  yjx2  + y2  > r0 


(4.2) 


where  uia  = —2.25,  r0  = 0.25  is  the  radius  of  the  circular  patches  (note  that  from  (4.1) 
A — + 1 when  <r  — > 0),  and  the  Gaussian  distribution  renders  the  vortex  compact  in  z.  We 
have  set  p to  give  an  initial  hydrostatic  balance  in  the  vertical  direction. 

The  vortex  goes  through  a brief  (on  the  scale  of  turn-around  time)  and  slight  adjust- 
ment and  then  evolves  into  an  equilibrium  shape  that  is  coherent  for  many  (we  have 
tested  up  to  120)  turn-around  times.  See  Figure  2. 


4.2.  Vortices  in  a Keplerian  shear 

Similar  to  the  case  without  shear,  we  start  with  an  initial  guess  that  is  close  enough  to 
the  equilibrium  state.  We  specify  our  initial  vorticity  function  as 


f wcexp  (--y)  ’ x2/a2  + y2/**2  ^ 1 

\ 0,  y/x2/a2  + y2/b2  > 1 


(4.3) 


where  u>0  = —0.625,  and  the  semi-axes  are  specified  as  a = 0.25  and  b = 1.0.  At  the 
midplane  2 = 0,  the  elliptical  vortex  patch  satisfied  exactly  the  Moore-Saffman  relation 
(4.1),  and  again  the  vorticity  strength  decreases  with  a Gaussian  distribution  away  from 
the  midplane.  Similarly  we  have  specified  the  initial  p such  that  the  flow  is  in  hydrostatic 
balance.  Note  that  we  have  specified  a lower  value  for  ui0.  This  is  because  J.A.B.  and 
P.S.M.  found  that  lower  vorticity  strength  (albeit  still  of  a moderate  Rossby  number) 
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Figure  3.  Elongated  three-dimensional  vortex  structure  in  a Keplerian  shear  (a  = 1.5).  Here  we 
show  the  isosurface  at  the  value  of  u>z  = —0.55,  and  at  the  time  t = 20  rori . Our  computational 
resolution  is  128  x 64  x 96  on  a computational  domain  of  8 x 2 X 6 (not  entirely  shown  in  this 
graph). 


x(H„) 


Figure  4.  Contour  plot  of  same  vortex  structure  as  shown  in  Figure  3 at  the  midplane  z = 0.  In 
this  plot  darker  color  designate  the  more  negative  values  of  wz , and  the  black  spot  at  the  center 
is  our  equilibrium  vortex  patch.  At  the  midplane  the  vortex  is  almost  exactly  Moore-Saffman. 

induces  fewer  thermodynamic  oscillations  in  the  vertical  direction  and  therefore  aids 
coherence. 

The  resulting  equilibrium  vortex  (after  a few  turn-around  times  from  the  initialization) 
is  shown  in  both  Figures  3 and  4.  In  Figure  3 we  show  an  isosurface  plot  of  the  vertical 
vorticity  ujz.  When  compared  with  Figure  2,  the  vortex  is  stretched  by  the  shear  and  has 
a elongated  shape  in  the  horizontal  directions  (Figure  4). 

Figure  5 illustrates  the  evolution  of  the  vertical  vorticity  profile  u>z  in  time.  The  vor- 
ticity quickly  adjusts,  on  turn-around-time  scales,  to  an  equilibrium  distribution  that  is 
maintained  for  a much  longer  period.  (We  have  tested  up  to  60  turn-around  times  in 
this  case.)  Note  that,  close  to  the  midplane,  the  vorticity  goes  through  small-amplitude 
oscillations.  This  is  because  at  z — 0 there  is  no  gravitation  and  the  flow  is  at  its  least 
stably  stratified. 

Finally  we  compare  our  results  (henceforth  denoted  as  “current”)  with  those  from  a 
sliding-box  implementation  (denoted  “sliding-box”)  that  has  been  developed  by  J.A.B. 
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Figure  5.  Time  evolution  of  vertical  vorticity  profile  u>z(z ) for  vortices  shown  in  Figure  3 and  4. 
The  vorticity  starts  with  a Gaussian  distribution  (solid),  goes  through  quick  adjustment  (dash) 
and  evolves  into  an  equilibrium  shape  (dot-dash  and  line-circle). 

and  P.S.M.  The  most  significant  difference  between  these  two  codes  lies  in  the  boundary 
conditions.  For  the  streamwise  x direction,  a natural  periodicity  is  assumed  for  both 
methods.  For  the  other  two  (transverse  and  vertical)  directions,  our  current  implemen- 
tation has  an  artificial  periodicity  in  z and  enforces  a no-flow  condition  in  y,  whereas  the 
sliding  box  formulation  assumes  a shear-convected  periodicity  in  y and  a no-flow  condi- 
tion in  z.  The  agreement  between  these  two  implementations  is  an  important  validation 
of  the  results  from  both. 

In  Figure  6 we  show  the  comparison  of  the  equilibrium  vorticity  ujz  along  the  three 
major  axes.  We  initialize  with  the  same  vorticity  distribution  that  is  given  by  equation 
(4.3).  The  velocity,  pressure  and  density  functions  are  obtained  satisfying  both  the  hy- 
drostatic balance  and  the  boundary  conditions  in  each  of  the  implementations.  The  flows 
then  evolve  and  we  compare  our  results  after  20  orbital  periods.  The  results  indicate 
excellent  agreement  between  the  methods. 

The  most  significant  difference  lies  in  the  vorticity  distribution  along  the  z axis.  This 
difference  can  be  better  illustrated  by  a contour  plot  of  ujz  along  the  x = 0 slice  (Figure 
7).  The  difference  in  wz(0,0,  z)  is  primarily  manifested  as  a smaller  vortex  size  in  the 
current  method.  Nonetheless,  the  equilibrium  shapes  of  these  vortices  agree  well  with 
each  other,  and,  most  strikingly,  the  black-and-white  stripe  patterns  emanating  away 
from  the  vortices  are  almost  exactly  reproduced  in  both  of  the  methods,  even  though 
the  boundary  conditions  are  very  different.  These  patterns  are  speculated  to  be  gravity 
waves.  Both  the  gravity  waves  and  the  difference  in  vortex  size  require  more  detailed 
study  in  the  future. 

5.  Future  Work 

In  the  present  paper  we  have  shown  the  numerical  implementation,  as  well  as  our  re- 
sults on  three-dimensional,  equilibrium  vortices.  We  have  compared  and  validated  our 
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Figure  6.  Comparison  of  equilibrium  vorticity  distribution  against  a sliding-box  calculation 
along  the  three  major  axes  of  the  vortex  structure. 


y(H,)  y(H„) 

Figure  7.  Comparison  of  equilibrium  vorticity  distribution  against  a sliding-box  calculation. 
Here  we  plot  the  function  w*(0 ,y,z),  and  the  darker  color  correspond  to  the  more  negative 
values.  The  very  black  patches  in  the  middle  are  our  equilibrium  vortex  structures.  Result  from 
the  current  method  is  shown  in  the  left  graph,  and  that  from  the  sliding-box  method  is  shown 
to  the  right. 


results  with  an  alternative,  sliding-box  method.  A few  of  the  issues  that  we  shall  be 
immediately  looking  into  next  are  as  follows.  First,  even  though  we  have  identified  a 
family  of  equilibrium  vortices  from  one  particular  form  of  initialization,  we  should  also 
seek  other  possible  forms  of  equilibrium  vortices.  For  example,  we  have  tried  the  combi- 
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nation  of  different  decaying  exponents  in  (4.3)  and  different  vorticity  strength  w0,  and 
found  equilibrium  vortices  under  certain  conditions.  (This  work  is  still  in  progress  and 
is  not  presented  here.)  Second,  we  have  preferred  a no-flow  condition  in  the  y direction, 
primarily  because  of  its  advantage  in  diagnosing  the  conservation  laws.  This  more  rigid 
boundary  condition  should  be  relaxed  or  modified  in  the  future  to  better  reflect  the 
physics  of  the  disk.  Finally,  the  purpose  of  the  current  numerical  implementation  has 
been  to  overcome  the  limitation  of  the  sliding-box  formulation,  and  eventually  simulate 
a radial  ( y ) thermodynamic  (potential-vorticity)  background.  Our  ultimate  goal  would 
be  to  generate  coherent  vortices  from  such  a (baroclinic)  background,  as  well  to  as  to 
study  the  dynamical  migration  and  deformation  of  such  vortices;  with  these  achieved  we 
will  be  ready  to  address  their  relevance  to  the  greater  problem  of  angular  momentum 
and  mass  transport  within  protoplanetary  disks. 
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